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ABSTRACT 

We present dynamical models based on a study of high-resolution long-slit spectra of the narrow- 
line region (NLR) in NGC 1068 obtained with the Space Telescope Imaging Spectrograph (STIS) 
aboard The Hubble Space Telescope (HST). The dynamical models consider the radiative force due 
to the active galactic nucleus (AGN), gravitational forces from the supermassive black hole (SMBH), 
nuclear stellar cluster, and galactic bulge, and a drag force due to the NLR clouds interacting with 
a hot ambient medium. The derived velocity profile of the NLR gas is compared to that obtained 
from our previous kinematic models of the NLR using a simple biconical geometry for the outflowing 
NLR clouds. The results show that the acceleration profile due to radiative line driving is too steep 
to fit the data and that gravitational forces along cannot slow the clouds down, but with drag forces 
included, the clouds can slow down to the systemic velocity over the range 100-400 pc, as observed. 
However, we are not able to match the gradual acceleration of the NLR clouds from '--^ to ~ 100 pc, 
indicating the need for additional dynamical studies. 

Subject headings: galaxies: kinematics and dynamics -galaxies: individual 

(NGC 1068) galaxies: Seyfert-AGN: [Oiii] emission lines -ultraviolet: galaxies 



1. INTRODUCTION 

Recent studies of the kinematics of the NLRs 
of Seyfert galaxies have taken advantage of the 
high resolution o f HST to map t he velocities of 
these r egions. (lEvans et al.l 119931: iMacchetto et al.l 



19941 iHutchings et al l 119981: iNelson et al l 120001: 



Crenshaw et al.l 120001: ICrenshaw &: Kraemeil l2000b 



Ruiz et al. ' 2001; ' Cecil et al.l l2002t iRuiz et al.l 12005 : 
Das et aL ll2005. 20061 and references therein). Although 



there have been a number of papers on the dynamical 
aspects of the NLR, most of these studies relied on 
ground-based data limited to spatial resolutions of > 50 
pc for even the most nearby AGN (see lKaiser et al.ll2000l 
and references therein). These studies relied primarily 
upon spatially integrated line profiles to understan d 
the dynamics of the NLR (|Schulzlll99fll : I Veilleuxl [l99l . 
but the problem arose that the emission-line profiles 
of the NLR can be explained by many different types 
of dynamical models, such as infall, rotatio n and 
outflows ( Capriotti et al. 198CLll981|; IVrtilekl 119831 : 
iKrolik fc Vrtilek,1984 : ,Vrtilek,19m 

With the launch of HST and its high angular resolu- 
tion {)". 1), the NLRs of Seyfert galaxies have received 
considerable attention. With the limited long-slit capa- 
bility of the faint object camera (FOG), and later the 
expanded capability of STIS, detailed constraints on the 
kinematics of the NLRs in Seyfert and other galaxies be- 
came possible. In turn, these kinematic studies provided 
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good diagnostics upon which dynamical analyses can be 
based. 

1.1. Previous Kinematics Studies of the NLR of Seyfert 

Galaxies 

The structure of the NLR resembles a bicone as ex- 
pected from a simple unified model of Seyfert galaxies, 
due t o collir nation by a th ick t orus (lAnton ucci fc Millen 
[1985[ ). Both lSchuld (|T990l ) and lEvans efa l. (1993) have 
modeled the NLR of NGC 4151 and found it to be consis- 
tent with a biconical geometry. In an HST study done on 
asample of Seyfert 1 and 2 galaxies, ISchmitt fc KinnevI 
(|1996l ) compared both Seyfert types to study their NLR 
morphologies, and found triangular structures in most 
of their Seyfert 2s, and circular structures in most of 
their Seyfert Is, consistent with the unif ied model and 
biconi cal s tructure for th e NLR (see also ISchmitt et all 
[2003alin ). IVeilleux et aP ((200l[l modeled the inner re- 
gions of the Seyfert galaxy NGC 2992, and found that 
the ionized gas can be fitted with a biconical structure. 

We have recently completed a study of the kin emat- 
ics of the NLRs in two Seyfert ga laxies ()Das et al.l 2005, 
hereafter Paper I, and iDas et al] 2006, hereafter Paper 
II) . Since our dynamical work in this paper is a direct ex- 
tension of the works of those two papers, we summarize 
our results below. In Paper I, kinematic models were de- 
veloped to match the emission-line velocities from high- 
resolution STIS spectra within ~ 400 pc of the central 
black hole of the Seyfert 1 galaxy NGC 4151. The NLR 
gas showed strong evidence of acceleration from ~ pc 
out to ^ 100 pc, then showed deceleration back to sys- 
temic velocity at ^ 400 pc with velocity roughly propor- 
tional to distance in each case. The maximum velocity 
of the outflowing gas at the turnover point (96 pc) was 
^ 800 km s~^ relative to the black hole. Based on our 
kinematic model, the NLR could be represented by a bi- 
cone, with inner and outer half opening angles of 15 and 
33° respectively, and inclination of ^ 45° with respect 
to the plane of the sky, consistent with previous kine- 
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matic work done on NGC 4 151 with different slit posi- 
tions (jCrenshaw et al.ll2000D . Some of the fainter NLR 
clouds showed evidence of backflow at the point where 
the clouds turnover in their velocities. The radio jet was 
found to have little effect on the kinematics of the NLR 
clouds, however there was some evidence of radial veloc- 
ity splitting of the fainter NLR clouds near bright knots 
in the radio jet. The brighter clouds were not accelerated 
by the jet. 

In Paper II, we developed a similar model for the 
Seyfert 2 galaxy NGC 1068 again with high-resolution 
spectra taken with the STIS aboard HST. With seven 
parallel slit positions covering the entire NLR, we ex- 
tracted radial velocity profiles and matched them with 
our newly developed 3-D biconical models. Our kinemat- 
ical models showed that the NLR gas accelerated out to 
~ 140 pc (the turnover point), then decelerated back to 
systemic velocity at a distance of ~ 400 pc from the cen- 
tral black hole. The maximum velocity of outflow of the 
gas was ^ 2000 km s~^, with respect to the black hole, 
and the model predicted and inner and outer half open- 
ing angles of the bicone of 20 and 40° respectively, and 
an inclination of 5° out of the plane of the sky. We used 
high resolution r adio maps of the NLR obtained from 
iGallimore et al.l (|2004f ) to search for jet/cloud interac- 
tions. Evidence showed that, similar to NGC 4151, the 
fainter NLR clouds were split near bright knots in the 
radio jets, whereas the brighter NLR clouds remained 
unaffected by the jet. 

Other Seyfert galaxies show similar flow patterns to 
NGC 4 151 and NGC 1068, suc h as the Seyfert 2 galaxy 
Mrk 3 (|Ruiz et al.ll200lL l2005t) . The radial velocities of 
these galaxies were matched with a common kinematic 
model with little variation in the parameter space, which 
begs the question, what are the physical processes involve 
that would cause such a similar flow pattern in both types 
of Seyfert galaxies? It is this question that motivated us 
to carry out this dynamical study. 

1.2. Previous Dynamic Studies of the NLR of Seyfert 

Galaxies 

In the flrst truly dynamical model of the NLR based 
on HST kinematics. lEverett fc Murray! ()2006[ ) attempted 
to fit the NLR velo cities in NGC 41 51, based on mea- 
surements done by iDas et all (|2005f) . They tested an 
isothermal Parke r wind model which assumes thermally 
expanding winds ()Parkeij|1965[ ). By assuming a spherical 
cloud geometry, they let the Parker wind drag along the 
embedded NLR clouds. As the Parker wind accelerates 
by thermal expansion and slowly loses heat by adiabatic 
cooling, the clouds are also accelerated to high speeds. 
They then let the Parker wind run into a low density 
ambient medium to slow the velocity of the wind, and 
hence the clouds, to the systemic velocity. The model 
explains the velocity profile of the NLR of NGC 4151, 
but suffers from the fact that an isothermal wind can- 
not be sustained out to large distances. Also, the mass 
profile of the SMBH plus galaxy, which determines the 
temperature profile for their model, is not exactly known 
for NGC 4151. Hence while their model was not success- 
ful on physical grounds, it is worth noting its relative 
success in matching the NLR kinematics. 

In this paper, the main question we want address is: 
how can we constrain the dynamics of the NLR in Seyfert 



galaxies with the detailed knowledge that we have gained 
from our kinematic studies? Here, we concentrate on the 
dynamics of the NLR in NGC 1068, since it has the best 
constraints. We will start with a simple construction of 
the enclosed mass function based on data from previous 
studies and eventually formulate a radiation pressure- 
gravity tug-of-war on the NLR clouds. The questions we 
will attempt to answer include the following: 1) If the 
NLR gas is in outflow, is radiation pressure really the 
best driving mechanism? 2) If the NLR gas is turning 
over its velocity and decelerating back to systemic, is 
gravity responsible for stopping the gas? 3) Can we fit 
the velocity profile of the data with a simple radiation- 
gravity law, or do we need to include another force (such 
as drag)? 

Our analysis applies to NGC 1068 in particular, but 
has relevance to all Seyferts in general that show signs 
of ga s outflow and subsequent deceleration (|Ruiz et al.l 
|2005| ). To test whether gravity is playing any role in 
stopping and turning back the gas velocity, we will con- 
struct a mass profile within ~ 10,000 pc from the nucleus 
of NGC 1068. With the mass profile in hand, we will test 
whether the gas kinematics are dominated by rotation. 
Such a test might prove fruitful for rotation in NGC 4151, 
which shows redshifts northeast and blueshifts southwest 
of the nucleus, but from a geometrical point of view, ro- 
tation of the gas will prove difficult to match the velocity 
field of NGC 1068, which shows blueshifts and redshifts 
on each side of the nucleus. Next we plot the escape ve- 
locity with distance to see if the gas should escape or not, 
given the velocities seen in the data, and whether or not 
the kinematics of the NLR can be dominated by gravita- 
tional infall. Then we concentrate on outflow assuming 
spherical symmetry and pure radial motion. First we de- 
termine if the deceleration of the gas can be attributed 
to the enclosed mass, regardless of the outward acceler- 
ating force. We then apply radiative line driving plus 
gravitational forces and compare the results to the ob- 
served velocity law of NGC 1068 derived from the kine- 
matic models. Finally, we introduce a drag force due 
to an ambient medium on the NLR clouds, in addition 
to radiation pressure and gravity, to determine if it can 
improve the flt to either the accelerating or decelerating 
portions of the observed velocity curve. 

2. BUILDING THE MASS PROFILE 

In building a model of the enclosed mass as a func- 
tion of distance from the central SMBH, we incorporate 
various subsystems into our mass proflle. These include 
contributions from the SMBH, the nuclear stellar cluster, 
and the bulge. For all these systems, we assume spheri- 
cal symmetry for simplicity. We have assumed that the 
stellar cluster and bulge extend all the way inward to the 
SMBH, which may overestimate the mass close in. The 
size of the stellar clu ster was estimated to be ^ 140 pc 
based on a study by iThatte et al.l (|1997l ) and the bulge 
was assumed to extend up to 10000 pc. Mass contribu- 
tion from the galactic disk of NGC 1068 was neglected 
because the galactic potential is dominated by the large 
bulge to at least 1000 pc, which is well beyond the extent 
of the NLR. 

2.1. The Black Hole Mass 
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NGC 1068 is one of only a few AGN that shows an 
edge-on disk of H2O maser emission close to the SMBH. 
The disk shows the signature of a rotational velocity 
curve, which can be used to deter mine the mass of the 
SMBH. 'GrccnhiU fc GwinnI (|1997l ) estimated the mass 
of NGC 1068 to be 1.5x10^ M© within 0.65 pc, based 
on the velocity field of the H2O maser emission observed 
with the VLBA and the VLA. We will therefore use this 
estimate for the mass of the SMBH. 

2.2. The Bulge Mass 

Since we are going out to ^ 10,000 pc, which encloses 
the NLR and the ENLR clouds, we need an accurate as- 
sessment of the total mass within this region. At small 
distances (< 1 pc), the SMBH is dominant, although in 
the case of NGC 1068, a concentrated stellar cluster is 
also providing substantial gravity close in. For an esti- 
mate of the bulge mass in NGC 1068, we rely on the work 
of Haring & Rix ( 2004). They found a tight correlation 
between black-hole mass and bulge mass for a sample of 
30 galaxies, including NGC 1068. They determined the 
bulge mass by modeling the bulge with the Jeans equa- 
tion in spherical form. They assumed the bulge to be 
isotropic and spherically symmetric, which might lead to 
an overestimation of the bulge mass; however they also 
neglect any contribution from dark matter, which would 
tend to underestimate the bulge mass. Their value for 
the bulge mass in NGC 1068 is 2.3x10^° Mq within a 
radius r = 3Re (3 effective radii), where Re = 3.1 ± 0.8 
kpc, a value taken from the sur face brightness deconvolu- 
tion of iMarconi fc HuntI (|2003D . The effective radius Re 
is defined to be such that half of the total light from the 
galaxy is predicted to be contained within the isoph otal 
ellipse that has area ttRI (jBinnev fc Merrifielcilll998D . 

Elliptical galaxies' and bulges' surface brightnesses can 
be well described by t he empirical formula developed by 
Ide Vaucoule"ini (|1948[) 



I{R) = lee 



-7.6692[(-jp-)4 -1] 



(1) 



where 1^ = I(-Re)- In the 1980s-1990s, a family of stel- 
lar density curves emerged that modeled both elliptical 
galaxies and bulges well. These curves are of the form 



p{r) = 



(3 - 7)M 



An 



r^{r + rf) 



4—7 ' 



(2) 



where 77 (in pc) is a scaling radius and M is the total mass 
of the bulge ()Dehnen„ 1993.1 . The parameter 7 determines 
different types of models, where the 7 = 2 c ases c orre- 
sponds to previous density m odels b y'Jaffe' (! l983[ ) and 
the 7=1 cases to models by iHernq uist (1990|). These 
models, when integrated over a spherical volume, yield 
the following enclosed mass profile: 



M(r) 



■iTTt^p{t)dt 



47r(3-7)Afr/ 
47r 



t^dt 



M 



3-7 



(3) 



For the special case when 7 = 1.5, the density pro- 
file of Equation [5] yields a surface density distribution 
that closely matches the de Vaucouleurs surface bright- 
ness profile of Equ ation [1] to wit hin 15% over nearly 4 
decades in radius ijDehnenI |1993[ ). Therefore we adopt 



the following form of the mass function as the bulge pro- 
file ^ ^ 

Mir) = M (^) . (4) 

We find a suitable value for ?; in Equation [5] by first 
defining the 'half-mass radius' ri: 



-M = M 
2 



(5) 



ri + 7] 

which yields the following relation for r/ 

7? = ri(2i-l). (6) 

iDehnenI ()1993[ ) has found a simple approximation for ^ 

that depends only slightly on 7. For a 7 < |, he found 
that 
R 

— « (0.7549- 0.004397+ 0.0032272-0.001827^). (7) 

ri 
2 

Therefore we can find a value for rj by using 7 — 1.5, Re 
given above, and Equation [SI We find that 77 ^ 2400 pc. 
We also know that M(3Re) = 2.3 xlO^" M©, so that 

n 1.5 



M = 2.3 X 10^7 



3 X 3.1 x 10^ 



3 X 3.1 X 103-1- 2400 
= 3.2 X 10^" Mq. (8) 
The bulge mass distribution can finally be written as 



Af (r) = 3.2 X 10 



10 



pc 



2400 



1.5 



M, 



0- 



(9) 



2.3. The Nuclear Stellar Cluster Mass 



It is known that NGC 1068 has a compact nuclear 
stellar cluster ~ 140 pc in radius (jThatte et al.l 119971 : 
ICrenshaw fc Kraemeil [2000al) . which contributes signif- 
icant mass to the total mass profile of the NLR of 
NGC 1068. Therefore this mass must be taken into 
account when deriving the mass profile. iThatte et al.l 
(1997) found a mass of 6.8x10^ Mq within 1" of the 
SMBH of NGC 1068, assuming a virialized, isotropic, 
and isothermal distribution of the stars. They used th e 
stellar velocity distribution (cr,) found in lDressleil (|1984f ). 
a value of 143 ± 5 km s^^ at ^ 1", to calculate the total 
dynamical mass within 1" of the nucleus of NGC 1068. 

2a^R 



G 



(10) 



where i? = 1" « 72 pc for NGC 1068. This mass includes 
contributions from the stellar cluster, the nucleus, and 
the bulge within a radius of 72 pc. Therefore to find just 
the mass from the stellar cluster Afsc, we took out the 
rest of the mass contribution from the total mass. 

Msc{72pc) = Mdyn{72pc) - M,„,bh - Mbulgei72pc) (11) 

Making the various substitutions we have 



Msc{12pc) = 6.8 X 10* 



3.2 X 10 



10 



1.5 X 10^ 
72 



72 + 2400 



1.5 



= 5.1 X lO'' Mp 



0- 



(12) 
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Based on the radial surface b rightn e ss profile presented 
in Fig ure 4 of iThatte et all (| 19971 ). iBeckert fc Duschll 
()2004f ) computed a power law consistent with the form 
S^,{r) oc to fit the data. They claim that if the 
profile traces the stellar mass distribution, then they 
would expect that the mass profile would have the form 
M*(r) oc r, given a spherical, isothermal distribution of 
stars in the cluster. Using their distribution we can es- 
timate the stellar mass function based on the condition 



kr, so that k — 



5.1x10° 



72 



Msc,72 = 5.1 X 10** M(T 

" Iq pc" 
therefore 

M,c(r) = 7.1 X lO^pc Mq. (13) 

Beyond the observed extent of the stellar cluster, we as- 
sumed that Msc{r) is a constant. Finally the total en- 
closed mass function for the NLR of NGC 1068 is given 

by Mtot{r) = M^mbh + Mtuige + Msc or 



Mtot{r) = 1.5 X 10^ + 7.1 X 10%c 



3.2 X 10' 



' pc 



2400 



1.5 



Mr, 



(14) 



3. DYNAMICS BASED ON GRAVITY 



A figure representing the total mass enclosed within rpc 
is shown in the top panel of Figure [1] The mass profiles 
for each contribution, the SMBH, bulge, and cluster are 
also shown in the figure. Close in toward the nucleus, the 
nuclear stellar cluster dominates up to its entire extent, 
while the bulge takes over from there. The black hole 
mass dominates at < 2 pc. The kink in the total mass 
curve at ~ 140 pc is because the stellar cluster was cut 
off abruptly at this location. We could have modeled 
the cluster to assume an exponential drop-off after 140 
pc, but this would not have contributed much to the 
gravitational force exerted by the total mass. 

3.1. Rotation 

With the total mass profile we calculate the circular 
rotational velocity and plot it as a function of distance 
as shown in the middle panel of Figure [TJ The rotation 
velocity only depends on the enclosed mass at radius r^c, 
and is given by the formula below 



(15) 



For demonstration purposes, we show the rotational ve- 
locities as if only each mass component was present, as 
well as the velocity for the total mass profile. The rota- 
tional velocity profile indicates to us that rotation cannot 
dominate the kinematics of the NLR of NGC 1068 be- 
cause many observed data points exhibit large velocities 
(> 1000 km s^'), whereas the rotation curve never ex- 
ceeds ^ 220 km s~^ at large distances (100 pc). Again 
the stellar cluster dominates the velocities up to ^ 300 
pc. 

3.2. Escape and Infall Velocity 

The escape velocity at a given distance r^c is calculated 
numerically from the formula 



Vir) = 



t'^ 



(16) 
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Fig. 1.— Top: The total enclosed mass profile of NGC 1068 
in solid black. The individual contributions to the total mass is 
also shown in this plot. They are the stellar cluster (dash-dot), 
the bulge (dash), and the SMBH (dot). The stellar cluster was 
modeled to have an extent of 140 pc. 

Middle: The total rotational velocity profile of the NLR of 
NGC 1068 in solid black color. The individual rotation compo- 
nents due to the various masses are also shown. 
Bottom: The escape velocity as a function of distance for the NLR 
of NGC 1068. 



Dynamics of the NLR in NGC 1068 



5 



based on the enclosed mass function, and is plotted in the 
bottom panel of Figure [TJ According to our kinematic 
model, the maximum velocity at the turnover radius rt = 
140 pc is 2000 km s^^ (Paper II); therefore Figure [1] tells 
us that the NLR clouds should have escaped after 140 pc, 
where the escape velocity is only ~ 500 km s^^. In the 
data however, clouds at 140 pc are at much higher ve- 
locity than escape velocity; yet after the turnover point, 
the clouds start to decrease their velocity and return to 
systemic. Thus, some force other than gravity is causing 
the clouds to decelerate at r > 140 pc. The infall ve- 
locity profile is equivalent to the escape velocity profile, 
except that the velocity vector is now directed inward. 
Therefore gravitational infall also cannot account for the 
faster moving clouds at 140 pc, and in general does not 
match the observed velocity profile (Paper II). 

3.3. Gravitational Drag 

To test the importance of the force of gravity alone 
on slowing down the outflowing NLR clouds, we give the 
clouds a maximum velocity at the turnover point and let 
gravity do the rest. In other words, we assume that there 
is no outward driving force after the turnover point and 
let the clouds coast under the force of gravity. The top 
panel of Figure [5] shows that with a maximum velocity 
> 1000 km s^^ at 140 pc, there is little deceleration with 
radius. However, with maximum velocities < 300 km s~^ 
gravity can slow the clouds down to rest, as seen in the 
bottom panel of Figure[2l The maximum velocities in the 
NLRs o f some Seyfert ga laxies are on the order of ~ 400 
km s~^ ()Ruiz et al.ll2005[ ). and gravitational deceleration 
may be important in these cases. The kinematic model 
of NGC 1068, however, shows maximum velocities of up 
to ^ 2000 km s~^, clearly out of the reach of gravita- 
tional deceleration. Gravity alone cannot slow down the 
outflowing clouds in this case and there must be some 
other force or forces involved. 

To compound the problem, suppose we let radiation or 
some other force push on the gas while gravity is trying 
to pull it back. In this case the gravitational deceleration 
will be even less. However, to further explore this issue, 
we generate a line-driven radiation model, with gravity 
competing to slow the gas down. This will eventually 
lead us in the direction of including a drag force to com- 
plete the analysis. 

4. RADIATION-GRAVITY FORMALISM 

The various mechanisms to push the gas out from close 
to the nucleus include a radiation pressure driven wind, 
a thermally driven wi nd, or a magn eto-hydrodynamic 
(MHD) wind (Crensh aw et al.l l2003f ). The latt e r two 
met hods are discussed i n mo re detail in lEverettl (j2005l ) 
and lEverett fc Murravl (|2006D . and were summarized in 
^1.2\ However, no dynamical model to date has led to a 
satisfactory description of the kinematics in NGC 1068 
or NGC 4151. The radiation driven wind mechanism 
is most efhcient when the momentum imparted to the 
gas is due to line-driving (bound-bound transitions), 
although bound-free and free-free electron transitions 
(Th omson scatterin g) also contribute to driving the gas 
out (jChelouche fc Nctzcr 2001). The radiation force is 
dependent on the ionization state of the gas, with lower 
ionization states more efficient due to the greater avail- 
ability of electrons in the bound states. If dust is mixed 
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10000 



Fig. 2. — The velocity slowing down with distance for an initial 
launch at 140 pc with 1000 km s~^ (top) and 400 km 8"^ (bottom) 
shown in solid lines for the clouds in NGC 1068. The other mass 
components' effect are shown in different linestyles represented by 
the key above the curves. The clouds barely slow down out to 
10,000 pc when launch with 1000 km s~^. However, gravity was 
able to slow down the clouds down when they were launched with 
300 km s-i. 



in with the NLR gas, then it will compete with the gas in 
absorbing ionizing photons and hence radiation pressure 
on the dust can become an important cont ributor to the 
veloc ities of the outflowing gas in the NLR (jDopita et al.l 
|2002| ). If the dust grains are electrically charged, they 
can drag the ionized gas along to similar velocities as the 
dust. In this section we ignore the effects of dust, which 
would only increase the radiative acceleration. Thus we 
consider a radiation driving mechanism coupled with the 
effects of gravity to find the velocity profile of the NLR 
gas. 

4.1. Building the Velocity Equation 

We start with the acceleration due to radiation on a 
point mass, 

LcjtM 

«W = -A — 5 ' (17) 
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where a is the acceleration, L is the bolometric lumi- 
nosity of NGC 1068, (JT is the Thomson scattering cross 
section for the electron, r is the distance, c is the speed 
of light, rUp is the mass of the proton, and M is the force 
multiplier. As mentioned above, to really drive the gas 
out efficiently, we need to incorporate other sources of 
opacity such as bound-bound and bound-free opacity in 
addition to those from Thomson scattering. These addi- 
tional opacities are included via the force multiplier, M. 
in Equation I17[ which is primarily a function of ioniza- 
tion parameter^ U for a gi ven spectral energy distribu- 
tion (jCrenshaw et al.ll2003l and references therein). 

The acceleration due to gravity per mass is simply 
given by 

a{r) = (18) 

where M{r) is the total enclosed mass within r parsecs, 
and G is the universal gravitational constant. Putting 
Equations [17] and [18] together we have 



a{r) 



LcttM GM{r) 



(19) 



Now we will have to rewrite the acceleration in terms of 
velocity as a function of radius, and then solve for v{r). 



(20) 



dv dr dv dv 

dt dt dr dr 

Therefore we can now write Equation [19] as a simple sep 
arable differential equation 



vdv 



LcttM , GM(r) , 
dr — dr. 



475-7-2 



(21) 



Substituting for the constants in Equation [21] and con- 
verting to appropriate units of km and pc, then 
integrating and setting the initial velocity to zero yields 
the following form for v{r) 



v{r) 



6840^44 



M 

t2 



8.6 X 10 



dt, (22) 



where L44 is luminosity in units of 10*^ ergs s~^ and 
M{r) is in units of Mq. The constraints on the lumi- 
nosity and the force multiplier are presented in the next 
section. 

4.2. Physical Constraints on the NLR of NGC 1068 

Since NGC 1068 is a Seyfert 2 g alaxy, we ca,nnot mea- 
sure its luminosity directly. From lPier et al.l (|1994D . the 
total luminosity of NGC 1068 is given by 



Lboi = 2.2 X 10^ 



0.01 ) \22MpcJ 



Lr. 



(23) 



where frefi is the fraction of nuclear flux observed as 
scattered radiation, D is the distance to NGC 1068, and 
Lq is the solar lumi nosity. The most uncertain term in 
Equation [23] is frefi- iPier et all (|1994l ) have summarized 
a range of values for frefi that have been determined pre- 
viously by several authors. The range in frefi spans a few 



The ionization parameter U is defined as U 



TOO ^ 



# of ionizing photons 



at the ionized face of the clouds, where hvQ 



# H atoms 

13.6 eV and nu is the hydrogen number density. 



Lhoi - 2.2 X 10^ 



= 2.4 X lO*'' ergs"\ 



orders of magnitude, from 0.0 01-0.05. They claim that 
the best estimate comes from iMiller et all (|1991[ ). who 
had determined a value for frefi to be 0.015, based on 
obser vations of fOml an d broad Hfi luminosity and their 
ratio. Pier et al] (|1994l ) concluded that frefi is probably 
within a factor of a few of 0.01, hence we have adopted 
a value for frefi of 0.015 because it is the " best" esti- 
mate a nd close to the average value adopted bv lPier et al.l 
(119941 ). We already know the distance to NGC 1068 as 
14.4 Mpc (jBland-Hawthorn et"allll997f) so the bolomet- 
ric luminosity of NGC 1068 is given by 

,11 ^o.oi5yY 14.4 

0.01 ) \22MpcJ ® 

(24) 

a value which could be uncertain by a factor of 0.3-15, 
depending on frefi- 

The emission lines arising from the NLR gas are 
best fitted with a two component photoionization 
model at e ach position, based o n HST/STIS long- 
slit spectra (jKraemer fc Crenshaw|[200Q) . Accord ing to 
iKraemer fc Crenshawl the ionization state of the gas 
ranges from U ~ 10~^'^-10~'^ ° for the two components 
but seems to vary little with distance. Using U and the 
SED for NGC 1068, we found the force multiplier to vary 
from Ai ~ 500-6000 for the front face (i onized face) of 
the c louds, based on CLOUDY models (| Per land et al.l 
[199^ . The mass function was derived in previous sec- 
tions, and L44 = 2.4, so we can now numerically solve 
Equation [22] for v{r), assuming A4 is constant with dis- 
tance. The results are presented in the next section. 

4.3. Radiation and Gravity Results 

Equation [22] has only two parameters that we can vary 
to find the velocity w as a function of distance r: the 
launch radius ri , and the force multiplier A4 . We plotted 
v(r) for various combinations of launch radii and force 
multipliers of the gas. The top panel of Figure [3] shows 
that with a launch radius of ri = 1 pc, the velocity of the 
gas increases with force multiplier, but quickly reaches a 
terminal velocity and does not slow down significantly, 
even out to ~ 10,000 pc. This is not a surprise, as previ- 
ous plots had shown that the mass is not enough to slow 
down the high-velocity clouds. In addition, the maxi- 
mum velocities for A4 > 500 are too high compared to 
the observations, indicating that we must increase the 
launch radius ri . 

The physical constraints on the NLR gas suggested 
that the force multiplier is larger than 500. However, we 
generated plots for smaller than 500 to see whether 
or not gravity can slow the clouds down, if the clouds are 
launched from 1 pc. In the middle panel of Figure [3] we 
see that with force multipliers < 40, gravity does slow 
the clouds down, but that the maximum velocities are 
too low to fit the data of NGC 1068. Such models may 
prove useful for Seyferts with lower outflow v elocities in 
their NLRs, such as N GC 4151 and Mkn 3 pas et all 
l200llR^fireraII[200^ 

We next fix the force multiplier to the lowest value of 
500 (because higher values of A4 will be even more prob- 
lematic to slow down the clouds) and vary the launch 
radius. These plots are presented in the bottom panel 
of Figure [3] The maximum velocity falls below the ob- 
served values at ri > 2 pc, and the final velocity does not 
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Fig. 3. — Top: The velocity profiles for several force multipliers 
A4 with a fixed launch radius ri = 1 pc. None of the curves show 
a significant decrease in velocity up to 10,000 pc. 
Middle: The velocity profiles for several lower force multipliers 
with ri = 1 pc. Curves with force multipliers less than 40 turn 
over before reaching 10,000 pc, although their maximum velocities 
do not fit the data of NGC 1068. 

Bottom: The velocity profiles for several different launch radii (top 
to bottom curves ri = 1, 2, 4, 8, 16, 32 pc) all modeled with a 
constant force multiplier of 500. The bottom curve was able to 
turn over and return to systemic at ~ 1000 pc. The maximum 
velocity and maximum distance reached for these curves do not fit 
the data. 



return to zero until we reach large launch radii. As the 
launch radius increases to ~ 32 pc, the maximum attain- 
able velocity of the outflowing clouds decreases, and the 
enclosed mass is finally able to turn the velocity around. 
However by that time, the maximum velocity is much 
lower than predicted by the kinematic model. Also, the 
launch radius and the maximum distance of outflow be- 
come too large to fit the data, since the NLR clouds 
generally launch from close in and drop in radial velocity 
within 400 pc (Paper II). Clearly these curves cannot fit 
the observed NLR velocities of NGC 1068. 

5. RADIATION, GRAVITY, AND DRAG FORCES 

The radiation-gravity interaction on the NLR clouds of 
NGC 1068 fails to reproduce its velocity profile. For rea- 
sonable parameters, the velocity increases quickly close 
to the nucleus and mostly remains constant over large 
distances regardless of launch radius. The maximum out- 
flow velocity is rather sensitive to launch radius and de- 
creases with increasing ri but the clouds' velocities never 
turn over and decrease. Fine-tuning ri and M outside 
of the range of reasonable parameters can lead to a de- 
celeration profile, but the resulting velocity profile and 
amplitude does not match the observed trend. Therefore 
we conclude that there must be additional forces at play 
in the NLR to account for the velocity profile that we see 
in the data. One such force that could explain the trend 
in the data is drag, whereby the clouds are slowing down 
in a more diffuse, hotter, and higher ionization medium 
(jCrenshaw fc Kraemeiir2000b[ ). 

5.1. Deceleration due to Drag 

The drag force exerted on a cloud by an ambient 
medium is 

^drag. cloud — Pmed{^cloud '^med) -^cloud: (^^) 

where Pmed is the mass density of the ambient medium, 
Vcioud and Vmed are the velocities of the cloud and 
medium respectively, a nd Adoud repre sents the cross sec- 
tional area of a cloud (|Everett fc Murrav 20061 Follow- 
ing Everett & Murray, we assume the clouds to be spher- 
ical with mass rricioud = ^'^RlioudPdoud, where Rdoud is 
the radius of a cloud, and pdoud is its mass density. The 
acceleration on the clouds due to drag is then 



^drag, cloud 



[Vcl oud 



'cloud 



z'^^cloudP cloud 



^p^cloud-^cloud 4 



, y^cloud 



(26) 



where n is the hydrogen number density and rrip is the 
mass of the proton. In our case, we assume that the 
velocity of the ambient medium is zero and the radius of 
the cloud remains constant. Therefore with Vmed ~ we 
can write 

l^med 3 2 (07\ 

O-drag, cloud — '^-jr; 'T'^cloud^ V^'J 

^^H, cloud 4 

where Nh, doud (« ndoudRdoud) is the hydrogen column 
density of the cloud, and the minus sign represents decel- 
eration. Combined with the radiative and gravitational 
acceleration of Equation[Tni the total acceleration on the 
clouds becomes 

LutM GM{r) n^ed 3 2 



O-tot 



A'Kr'^c 



Nh,cIo 



"clouds 



(28) 
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and in differential form is 

LgtM , GM{r) , 
vdv — - — 5 dr — dr 



v^dr. (29) 



Substituting the various constants and converting to 
units of pc and km s^^, we have the following insep- 
arable differential equation to solve 



dr 



rir. 



4.3 X 10~ 



-2.3xf0 -ir^v{r), 
where N20 = iV/^/lO^^cm-^. 



(30) 



5.2. Constraints on Cloud and Medium Densities for 
NGC 1068 

Both NGC 1068 and NGC 4151 show evidence for 
highly ionized gas extended throughout t heir NLRs, 
based on Chandra X-ray Observatory images (jOgle et al.l 
HOOO, 2003). Thus, we assume an ambient medium that is 
highly ionized, with ionization parameter U « 10. The 
ionization parameter, which is inversely related to the 
density and radius, can be written as follows 



U oc 



r^riH 



(31) 



Therefore if the ambient and the NLR clouds see the 
same ionizing luminosity (Lion) at a particular distance 
r from the source, we can write 



rf'med — ri^iQini 



med 



(32) 



iKraemer fc Crenshaw! ()2000[ ) provided good constraints 
on the parameters on the right side of Equation [321 The 
densities of the NLR clouds are almost constant out to 
large distances from the nucleus with a typical value of 
ricioud ~ 10^ cm~'^ for clouds with Udoud ~ 10"'^. If we 
substitute these estimates in Equation 1321 we will have 
a typical estimate for Umed as follows 



n^med 



10^ 



10 



-3 



10 



1 cm ^. 



(33) 



IKraemer fc Crenshaw! found column densities for the 
NLR clouds in the range Nh — 10^^-10^^ cm~^, which 
corresponds to N20 = 0.1-10. Since we are interested 
in the ratio "j^^^ , varying either parameter while keep- 
ing the other constant will result in the same curves. In 
this paper, we choose to keep iV2o constant at the aver- 
age value of 1 and vary Umed, because Umed is the most 
unknown quantity. We already know that A4 can take 
values from 500-6000, so Equation !50! can now be solved 
numerically for v{r) with various values of the launch 
radius (ri), force multiplier {Ai), and the ratio ^j^- 
We used Mathematica v5.2, which employs the most ef- 
ficient choices among various flavors of Runge-Kutta al- 
gorithms, to solve for v(r). The results are presented in 
the next section. 

5.3. Radiation, Gravity and Drag Results 

The top panel of Figure |3] represents several plots with 
a force multiplier of 500, a launch radius of 1 pc, a col- 
umn density of iV2o = 1, and various densities of the 
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Fig. 4. — Top: The velocity profiles for force multiplier 500, a 
launch radius of 1 pc, and varying medium densities n^ed shown 
by the numbers. 

Bottom: The velocity profiles for a force multiplier of 6000, a 
launch radius of 20 pc, and varying medium densities similar to 
the top panel. 

ambient medium. The figure shows that when launched 
at 1 pc, the gas accelerates to a maximum velocity in- 
side 10 pc, then slows down again. The velocity slows 
down faster with increasing drag forces, as measured by 
increasing densities of the ambient medium. However the 
point of maximum velocity is too close in to match the 
data, which has a turnover velocity at ^ 140 pc. The 
curve with medium density n^ed = 0-33 shows gradual 
deceleration out to ^ 500 pc, but its maximum velocity is 
slightly too high. Any other curve has either too high ve- 
locity at turnover or the velocity drops too quickly. Note 
that the velocity of the data we are trying to match has 
- 2000 km s-i at a turnover of - 140 pc for NGC 1068. 
With force multipliers higher than 500, the same trend 
as in Figure [4] is seen in the velocity model, except that 
the velocities are much higher (> 6000 km s^^). Again, 
the curve that seems best to match the deceleration part 
in the data is with medium density nmed = 0.333. The 
rest of the curves decelerate too quickly or too slowly. 
The major problem with all of these curves is that the 
velocity turnover point is much closer to the nucleus than 
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TABLE 1 

Parameters used to generate the drag model shown in 
Figure [5] 



^max 


^iniicr 


^outcr 




(pc) 


(deg) 


(deg) 


(deg) (deg) 


450 


10 


40 


5 57.8 Equation l30l 



Zmax is the distance from the center to one end of the bicone, 
along the z-axis. 

Sinner is the inner opening angle of the bicone. 

Souter is the outer opening angle of the bicone. 

iaxis is the inclination of the bicone axis out of the plane of 

the sky, with positive inclination implied by the north bicone 

closer toward the observer. 

Praxis is the position angle of the bicone axis in the plane of 
the sky. 

^lain is the velocity law used in mapping the velocity field onto 
the bicone. 

140 pc. 

Another way to decrease the overall velocity is to in- 
crease the launch radius, and tweak the force multiplier 
and ambient density. We tried a launch radius of 10 pc 
with force multiplier ranging from 500-6000. The maxi- 
mum velocity drops as expected, reaching ^ 700 km s^^ 
for a force multiplier of 500 and climbs to ~ 2500 km s^^ 
when we increase the force multiplier to 6000. We no- 
ticed that the turning point increases by a factor of 5 as 
we increase the launch radius from 1 to 10 pc. However 
the turnover point is still too low to fit the data well. In 
order to have a good fit, we first need to maximize the 
launch radius and then tweak the force multiplier and 
medium density to get the closest match to the data as 
possible. 

The resolution in our data is ~ 10 pc and we should 
not launch much beyond this distance since we see clouds 
close to the SMBH with near zero velocity. Therefore, in 
the bottom panel of Figure |4l we present a plot with a 
launch radius of ri — 20 pc, a force multiplier of 6000, 
and a column density of N20 = 1 . The maximum velocity 
reached is ~ 1700 km with Umed = 10. The best 
curve to represent the data seems to be the one with 
nmed = 0.33, whose maximum velocity is 1500 km s~^, 
although the turnover point at ^ 60 pc is still too low 
according to our kinematic model. Furthermore, we have 
had to tweak the launch radius and force multiplier to 
very specific values, such that only a narrow range of the 
observed values gives a reasonably decent fit. However, 
to directly test this velocity profile, we generate biconical 
models similar to our previous kinematic models (Paper 
II) , to determine whether a more dynamical velocity law, 
rather than the simple linear law, can reasonably match 
the data. 

5.4. Dynamical Fit 

We applied our kinematic models of Paper II with our 
dynamical velocity law to the data. Previously we had 
used a simple kinematic velocity law, which is based on 
the relation v = kr. Now instead we substitute the ve- 
locity relation based on Equation 1301 which represented 
a more physical situation. A model with the new veloc- 
ity law of Equation 1301 is presented in Figure EJ using the 
input parameters from Table [U and using the best-fit 
curve from the bottom panel of Figure 21 with A^20 = 1 
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Fig. 5. — The model of slit 4 generated with input parameters 
from Table [T] and the best-fit velocity profile from the bottom 
panel of Figure |4] with n^g^^ = 0.33. Clearly this model is a poor 
match to the data. 

with an inner and outer half opening angle of 10° and 40° 
respectively, and inclination of 5°, and position angle of 
57.8°, and a half-size of 450 pc. Points on the bicone were 
assigned velocities according to Equation [3D1 The cen- 
ter slit (slit 4, see Paper II) is used here for comparison. 
We extracted the velocities from the model at a position 
equivalent to slit 4 and those are shown in the shaded 
regions in Figure [5] The data from slit 4 are shown in 
small triangles. 

The model represents a poor fit to the data from slit 4 
(the center slit) of NGC 1068. This was expected from 
looking at the previous velocity plots, as the turnover 
point was too low, the launch radius was a bit large, 
and the velocity profile did not resemble our kinemati- 
cally derived linear profiles. We have varied the input 
parameters from Table [l] but this makes very little im- 
provement. The turnover point, starting distance, and 
maximum velocity of the bicone cannot be varied with- 
out changing the drag parameters, as these are implic- 
itly defined in Equation [301 The thickness of the bicone 
can be increased to accommodate more data points, but 
the bicone will show lots of unnecessary shaded regions. 
Changing the maximum extent of the bicone Zmax will 
have absolutely no effect on the shaded region, except in- 
terrupting the shaded regions before or continuing them 
beyond ± 6". That leaves us with only two parameters 
to vary, the inclination of the bicone axis, and its posi- 
tion angle in the sky, and varying these two alone did 
not fix the model. 

6. CONCLUSIONS 

With radiation pressure driving the NLR clouds, their 
velocities will accelerate very quickly, within a few par- 
sees of the nucleus, assuming the clouds are indeed 
launched close in. With the introduction of the drag 
forces, the overall velocities are lowered, but even so the 
velocities reach maximum too quickly. The data suggest 
that the clouds are gradually increasing their velocities 
to a maximum at about 140 pc. This gradient in the ve- 
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locity cannot be simply accounted for by radiative forces 
driving the clouds. It seems therefore, that radiation 
pressure may not be the only driving mechanism for the 
NLR clouds, or that other forces are involved to steer the 
clouds to that particular gradient. The velocity profile 
shown in the data resembles one with a linear or 'Hubble 
flow' law. The inclusion of a drag force only serves to re- 
duce the overall high velocities due to radiation driving, 
but the maximum velocities are reached very close to the 
nucleus, too close to match the data effectively. 

Gravitational forces alone cannot stop the fast mov- 
ing clouds observed in the NLR of NGC 1068. With 
velocities as high as 1500 km s~^ at 140 pc, the clouds 
should have escaped the NLR. To compound the prob- 
lem, when radiation forces are added, the cloud veloci- 
ties are boosted even more. Yet we see clouds that are 
slowing down and gradually reaching systemic velocity. 
Therefore, the data suggest that there is a powerful force 
dragging on the clouds to slow their velocity. 

The drag force that we introduce can have a significant 
effect on the clouds' velocities. We can conclude there- 
fore that the drag forces are a strong competitor to the 
radiative forces, strong enough to bring the clouds to a 



halt even close to the nucleus, depending on the column 
densities of the outflowing clouds and the densities of the 
ambient medium. However the overall velocity profiles 
generated with radiative, gravitational, and drag forces 
do not match the data for NGC 1068. Assuming that 
the mass profile of NGG 4151 is similar to the one for 
NGC 1068, it will prove difficult to match its observed 
velocity profile, because the same linear trends are seen 
in the velocity of t he outflowing clo uds. The same can 
be said for Mkn 3 (jRuiz et al.ll2001[ ). 
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